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Chapter  2 

Frequency  Domain  Wave  Models  in  the  Nearshore  and  Surf  Zones 

James  M.  Kaihatu 

Ocean  Dynamics  and  Prediction  Branch,  Oceanography  Division  (Code  7322) 

Naval  Research  Laboratoiy,  Stennis  Space  Center,  MS  39529-5004 

1.  INTRODUCTION 

In  deep  water  (kh>>  U  where  k  is  the  wave  number  and  h  the  water  depth),  second-order  wave 
nonlinearity  can  be  described  as  a  small  correction  to  the  underlying  linear  wave.  Perturbation  ex¬ 
pansions  in  wave  steepness  €  =  ka,  where  a  is  the  wave  amplitude,  are  used  (Phillips,  1960),  and  at 
second-order  only  non-resonant  (bound)  waves  are  possible  among  triads  of  wave  frequencies.  Thus 
the  interacting  waves  with  the  frequency-vector  wave  number  combination  (cox, ki)  and  {(02,k2) 
excite  secondary  waves  at  {o}\  +  a)2,  ki  -f  k2),  but  these  secondary  wave  amplitudes  always  remain 
small  relative  to  the  primary  amplitudes.  At  the  next  order,  resonant  interaction  occurs  between  quar¬ 
tets  of  waves,  with  the  resultant  slow  energy  exchange  between  the  interacting  waves. 

In  shallow  water  {kh  «\)  waves  become  less  dispersive  and  more  collinear,  and  triads  of  waves 
at  second-order  begin  to  more  closely  satisfy  the  resonant  conditions  for  wave  interaction.  The  per¬ 
turbation  solutions  of  finite  depth  do  not  apply  in  the  nearshore,  since  significant  energy  transfer 
occurs  over  much  shorter  distances  (0(10)  wavelengths)  than  in  deep  water.  The  Ursell  number 
Ur  =  afk^h}  (Ursell,  1953)  is  the  typical  measure  for  the  validity  of  these  perturbation  solutions, 
which  are  only  applicable  if  Or  <<  1.  Though  the  resonant  conditions  between  triads  are  only  ex¬ 
actly  satisfied  in  the  collinear,  non-dispcrsive  limit,  the  nonlinearity  inherent  in  shoaling  waves  in  the 
nearshore  is  strong  enough  for  significant  energy  transfer  to  occur  at  near-resonance  (Bryant,  1973). 
Recourse  is  often  made  to  the  Boussinesq  equations  (Peregrine,  1967)  for  simulation  of  nonlinear 
energy  transfer  in  shallow  water,  as  they  are  valid  for  Ur  =  0(1),  where  weak  nonlinearity  and  weak 
dispersion  are  balanced. 

1.1.  The  Frequency  Domain 

One  model  framework  which  has  been  used  in  simulating  ocean  wave  propagation  in  the  nearshore 
has  involved  the  application  of  Fourier  transforms  to  the  dynamical  equations  governing  the  propa¬ 
gation.  This  transformation  involves  imposing  the  following  constraint  on  the  dependent  variable  of 
these  equations  (usually  the  free  surface  r;) 

N 

«=1 

where  o)n  is  the  /zth  frequency  in  the  spectrum,  is  the  total  number  of  frequency  components  in 
the  spectrum,  is  a  complex  Fourier  amplitude  and  c.c.  denotes  complex  conjugate.  Assumption 
of  temporal  periodicity  is  a  natural  application  to  ocean  waves. 

The  frequency  domain  format  allows  explicit  detail  of  nonlinear  wave-wave  interaction  and  wave 
transformation  properties.  Nonlinearities  in  the  equations  appear  as  products  of  amplitudes  at  discrete 
frequencies  in  the  spectrum,  which  can  then  be  investigated  in  detail.  Since  the  resulting  equations  are 
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in  tenns  of  evolving  amplitudes  rather  than  the  free  surface,  spatial  resolution  requirements  are  usu¬ 
ally  less  restrictive  than  in  time  domain  models.  Overall  computational  time,  however,  is  a  function 
of  the  number  of  frequency  components  kept  in  the  simulation,  whereas  (outside  of  ensuring  suf¬ 
ficient  resolution  for  the  shortest  waves)  this  is  not  germane  for  time  domain  models.  Additionally, 
there  is  often  a  disconnect  between  properties  of  a  time  domain  model  and  those  of  the  corresponding 
frequency  domain  models,  A  good  example  is  seen  in  Rygg  (1988),  who  used  a  time  domain  model 
of  the  classical  (shallow  water)  Boussinesq  equations  of  Peregrine  (1967)  to  simulate  intermediate 
depth  cases  of  laboratory  wave  propagation  experiments  successfully.  Similar  experiments  with  cor¬ 
responding  frequency  domain  models  (Liu  et  al.,  1985,  as  used  by  Kaihatu  and  Kirby,  1995)  have 
proven  less  favorable. 


2.  CLASSICAL  BOUSSINESQ  MODELS  IN  THE  FREQUENCY  DOMAIN 


The  Boussinesq  equations  can  be  derived  from  either  the  Euler  equations  (Peregrine,  1967)  or  the 
boundary  value  problem  for  water  waves  (Mei,  1983).  In  shallow  water,  it  is  reasonable  to  assume 
that  vertical  velocities  in  the  water  column  are  much  smaller  than  horizontal  velocities.  This  imposes 
the  following  scales  on  the  independent  variables 


(x,y)  , 

L  ^  ^ 


(2) 


where  L  is  a  characteristic  wavelength,  ho  a  characteristic  water  depth,  and  the  primes  denote  di¬ 
mensionless  variables.  These  scales  are  then  applied  to  the  physical  quantities 


(«'.  u')  = 


ho 


t«'  = 


oL^/  gho 


(3) 


(4) 


where  a  is  a  characteristic  amplitude,  (m,  v,  w)  are  the  water  particle  velocity  components,  p  is 
the  pressure  and  p  is  the  fluid  density.  When  substituted  into  the  Euler  equations,  the  following 
dimensionless  parameters  become  evident 


S=^  (5) 

ho 

which  are  measures  of  frequency  dispersion  and  nonlinearity,  respectively.  The  Boussinesq  equations 
can  be  derived  by  assuming 

0(M^)«0(5)«C»(1)  (6) 


Using  the  scaled  Euler  equations,  Peregrine  (1967)  derived  the  Boussinesq  equations  for  a  varying 
bathymetry 

tit  +  Vih  +  r))u=0(ti*,Sp?',S^)  (7) 

h  h^ 

u,  +  n  •  Vu  +  gVf)  =  -  V  [V  •  (hu,)]  -  —V  [V  •  u]  +  Sp?,  S^)  (8) 

where  u  is  the  depth-averaged  velocity  vector.  The  quadratic  nonlinear  terms  in  the  equation  above 
represent  the  lowest-order  nonlinearity  of  0(8).  Application  of  Fourier  series  to  these  terms  requires 
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special  treatment  (see  Mei,  1983),  and  thus  gives  rise  to  the  triadic  cross-spectral  energy  transfer 
which  is  the  manifestation  of  nonlinearity  in  the  frequency  domain.  Freilich  and  Guza  (1984)  de¬ 
rived  frequency  domain  models  from  the  one-dimensional  form  of  these  equations.  The  first  (the 
“consistent”  model)  can  be  written 


in^k^h^  ^ 
- 2 - 


+ 


3ink 

W 


N-n  > 

AtA„-i  +  2  ^  AfA„+i 
1=1  ) 


(9) 


where  An  are  complex  amplitudes  of  the  free  surface  and  asterisks  denote  complex  conjugate. 
Freilich  and  Guza  (1984)  solved  the  equation  in  terms  of  coupled  amplitude  and  phase  equations 
rather  than  the  complex  amplitudes  seen  in  equation  (9).  The  second  model  (the  “dispersive”  model) 
was  also  derived  from  the  standard  Boussinesq  equations,  but  does  not  contain  the  phase-shifting 
dispersive  term  (third  term  in  equation  (9)).  Instead,  the  weak  dispersion  is  incorporated  through  the 
use  of  the  dispersion  relation  for  the  Boussinesq  equations 


(10) 


One  consequence  of  the  use  of  this  dispersion  relation  is  that  the  wave  number  kn  is  no  longer  a  linear 
function  of  co.  Thus,  the  interacting  amplitudes  (An,  An±/,  Aip/),  while  resonant  in  frequency,  are  in 
near  resonance  in  wave  number.  Freilich  and  Guza  (1984)  then  compared  both  models  to  field  data, 
using  offshore  wave  spectra  to  initialize  the  model  and  ably  demonstrating  the  utility  of  frequency 
domain  models  to  nearshore  wave  propagation  problems.  Their  comparisons  of  wave  spectra  showed 
that  the  dispersive  shoaling  model  performed  slightly  better  than  the  consistent  model;  however,  both 
models  clearly  deviated  from  the  data  in  the  higher  frequency  range,  where  kh  was  no  longer  small. 

Two-dimensional  frequency  domain  models  of  both  the  Boussinesq  equations  (7)  and  (8)  and  the 
Kadomtsev-Petviashvili  (KP)  equations  (Kadomtsev  and  Petviashvili,  1970)  were  developed  by  Liu 
et  al.  (1985)  in  the  form  of  parabolic  models,  which  are  formulated  based  on  the  assumption  that  the 
angle  between  the  wave  direction  and  the  jc-axis  of  the  grid  is  small.  Kirby  (1990)  developed  angular 
spectrum  models  based  on  the  Boussinesq  equations  of  Peregrine  (1967).  Periodicity  in  both  time 
and  longshore  direction  was  assumed,  thus  imposing  resonant  interaction  among  longshore  wave 
number  modes  as  well  as  frequency  modes. 


3*  EXTENDED  BOUSSINESQ  MODELS  IN  THE  FREQUENCY  DOMAIN 

One  fundamental  problem  with  frequency  domain  models  of  the  classical  Boussinesq  equations 
is  their  lack  of  applicability  in  deeper  water  than  that  for  which  the  shallow  water  theory  is  valid. 
Recent  efforts,  beginning  with  Witting  (1984),  have  focused  on  improving  the  deep  water  behavior 
of  Boussinesq  models  such  that  their  linear  properties  (dispersion,  shoaling,  etc.)  better  mimic  those 
of  fully-dispersive  linear  theory  for  a  wide  range  of  water  depths.  McCowan  and  Blackman  (1989), 
Madsen  et  al.  (1991)  and  Nwogu  (1993)  represent  some  of  the  first  attempts  to  improve  time  domain 
Boussinesq  models  in  this  regard;  the  resulting  models  were  generally  dubbed  “extended”  Boussi¬ 
nesq  models  because  their  linear  properties  were  extendable  to  intermediate  and  deep  water.  Madsen 
et  al,  (1991)  added  terms  to  the  classical  Boussinesq  momentum  equation,  multiplied  by  a  free  pa¬ 
rameter,  which  would  be  zero  in  shallow  water  but  have  significant  effect  in  deeper  water.  This  free 
parameter  was  tuned  via  Pad6  approximations  so  that  the  dispersion  relation  of  the  equations  would 
compare  favorably  to  that  of  linear  theory  for  a  wide  range  of  depths.  Madsen  and  S0rensen  (1992) 
extended  the  Madsen  et  al.  (1991)  model  to  include  varying  bathymetry.  Madsen  and  S0rensen  (1993) 
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investigated  frequency  domain  formulations  of  the  model  of  Madsen  et  al.  (1991)  for  wave  evolution 
over  a  flat  bottom,  and  sloping-bottom  extensions  of  this  equation  became  the  basis  for  further  de¬ 
velopment  in  the  frequency  domain  (Eldeberky  and  Battjes,  1996;  Kofoed-Hansen  and  Rasmussen, 
1998;  Becq-Girard  et  al.,  1999).  The  equations  of  Madsen  et  al.  (1991),  and  their  various  nonlin¬ 
ear  and  dispersive  enhancements,  have  been  analyzed  extensively  by  Schaffer  and  Madsen  (1995), 
Madsen  and  Schaffer  (1998)  and  Madsen  and  Schaffer  (1999). 

In  contrast,  but  to  the  same  end,  Nwogu  (1993)  used  the  velocity  variable  at  an  arbitrary  location  in 
the  water  column  (rather  than  the  depth-averaged  velocity  as  in  the  classical  Boussinesq  equations)  as 
a  basis  for  deriving  extended  Boussinesq  equations  from  the  inviscid  Euler  equations.  The  resulting 
equations  contained  higher-order  terms  in  both  the  continuity  and  momentum  equations,  and  are 


»?/  *f  V  '  [(/i  +  77)ua]  +  V  - 


hV[V  ■  (hUa)] 


n„,  +gVn  +  (Ua  ■  V)U<,  +  Z„  j  ^  V(V .  u<,f)  +  V[V  •  (/iu„,)]  j  = 


(11) 

(12) 


where  is  the  horizontal  velocity  at  a  location  Za  in  the  water  column.  The  dispersion  relation  of 
this  set  of  equations  is  found  by  isolating  the  linear  terms  and  substituting  in  a  periodic,  progressive 
wave,  leading  to 


1  --ot{kh)^ 


(13) 


where  a  is  a  free  parameter  related  to  Za  by 


This  free  parameter  oc  is  then  best-fit  to  that  of  fiilly-dispersive  linear  theory  for  a  wide  range  of 
water  depths.  Nwogu  (1993)  determined  that  a  =  —0.390  was  the  best-fit  parameter  value  for  the 
range  0  <h/Lo  <  0.5,  where  Lo  is  the  deep  water  wavelength.  This  value  of  ot  corresponds  to 
Za  =  -0.522A. 


3.1.  Frequency  Domain  IVansformation  of  the  Equations  of  Nwogu  (1993):  Linear  Properties 

Usually  the  first  step  undertaken  in  a  frequency  domain  transformation  is  to  combine  the  con¬ 
tinuity  and  momentum  equations  into  one  via  the  use  of  first-order  substitutions.  Noting  the  extra 
dispersive  terms  in  both  the  continuity  equation  (11)  and  momentum  equation  (12),  Chen  and  Liu 
(1995)  commented  on  the  difficulty  in  determining  a  frequency  domain  form  of  the  equations  such 
that  the  linear  dispersion  relation  (see  equation  (13))  would  remain  applicable  to  the  resulting  equa¬ 
tion.  Later,  Kaihatu  and  Kirby  (1998)  determined  a  series  of  first-order  substitutions  which  would 
lead  to  a  set  of  equations  retaining  the  original  dispersion  relation. 

To  illustrate  the  difficulty,  we  reduce  equations  (11)  and  (12)  to  their  linear,  one-dimensional  form 
for  a  flat  bottom 


r]t-^huax  + 


h^Uaxxx = 


0 


(15) 


Uat  +  grjx  +  oth^Uxxt  =  0 


(16) 
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We  follow  the  procedure  of  Liu  et  al.  (1985)  to  formulate  the  frequency  domain  model.  The  first  step 
involves  combining  the  continuity  (15)  and  momentum  (16)  equations.  We  make  use  of  the  following 
first-order  relations 


rit  =  -huax 


(17) 


Uat  -  -gVx 


(18) 


We  then  take  the  time  derivative  of  equation  (15),  the  jc -derivative  of  equation  (16),  and  combine  the 
resulting  equations.  We  then  use  equation  (18)  to  eliminate  Ua  in  favor  of  rj.  This  results  in 

ntt  -  ghrixx  +  gh^nxxxx  “  h^rixxxx  =  0  (19) 


To  obtain  the  linear  dispersion  relation,  we  substitute 

t)  =  (20) 

into  equation  (19)  and  obtain 

=  ghk^  1^1  -  (21) 

which  is  essentially  the  linear  dispersion  relation  to  weakly-dispersive  Boussinesq  theory  to  within 
a  binomial  expansion.  The  substitution  sequence  used  to  collapse  the  two  equations  did  not  retain 
the  dispersion  relation  of  the  original  equation.  Schaffer  and  Madsen  (1995)  addressed  this  issue  by 
applying  differential  operators  of  O(fjtP')  (multiplied  by  free  parameters)  to  the  equations  of  Nwogu 
(1993)  and  then  used  a  Pad6  [4, 4]  expansion  to  determine  the  set  of  parameters  which  best  fit  the 
linear  dispersion  and  shoaling  characteristics  from  linear  theory,  with  the  results  of  Nwogu  (1993) 
representing  a  subset  of  the  parameters.  In  contrast,  Kaihatu  and  Kirby  (1998)  used  a  different  series 
of  substitutions  to  retain  the  dispersion  properties  of  the  original  equation;  this  is  examined  here.  If 
we  had  taken  the  time  derivative  of  equation  (15),  and  then  used  equation  (18)  to  replace  Ua  with  rj, 
we  would  have  obtained 


mt+huaxt-g(ci  +  j^h^Tixxxx=0  (22) 

We  then  multiply  equation  (16)  by  ^  and  substitute  the  time  derivative  of  equation  (17)  to  elinainate 
Ua .  Substituting  the  result  into  equation  (22)  results  in 


mt  -ghrjxx-oth^rixxtt  -  gh^ 


fjxxxx  —  0 


(23) 


It  can  be  shown  that  the  linear  dispersion  relation  of  equation  (23)  is  equation  (13),  the  original 
dispersion  relation  of  Nwogu  (1993). 

The  complicated  substitution  sequence  required  to  retain  the  linear  dispersion  relation  also  affects 
the  shoaling  behavior  of  the  frequency  domain  model.  To  examine  this,  we  return  to  the  derivation 
of  equation  (23),  but  retain  bottom  slope  terms.  Performing  the  same  series  of  substitutions  and 
neglecting  hxx  and  (hx)^  terms  leads  to 


Yitt  ““  g{hrix)x  +  2ahhxr}xtt  +  C£h^r}xxtt  —  gh^(5a  -h  2)kxr}xxx  —  gh^ 


JJxxxx  =  0 


(24) 
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Substituting  equation  (20)  into  equation  (24)  leads  to 

Ax  +  WA^O 

(25) 

where 

Ekx  4-  Fhx 

G 

(26) 

E  =  gh  +  oP’h^a  -  6gh^  ^a  + 

(27) 

F  =  gk  +  Tacohh  -  ghHSa  +  2)k^ 

(28) 

G  =  2  ^kh  +  oP'h^ak  -  2gk^h^  ^a  + 1  j  j 

(29) 

Though  the  derivation  appears  to  be  fairly  straightforward,  it  will  be  shown  that  the  linear  shoaling 
term  (equation  (26))  compares  very  poorly  to  that  of  linear  theory.  Further  analysis  reveals  that  the 
balance  between  the  rjttx  and  g(hTjxx)x  terms  governs  the  effectiveness  of  the  wave  shoaling  relation. 
Kaihatu  and  Kirby  (1998)  addressed  this  by  adding  the  following  term  to  the  equation 

-  ghrjxx)x  =0  (30) 

which  is  true  at  lowest  order.  In  this  equation  fi  is  r  free  parameter  to  be  optimized.  This  changes 
equation  (24)  to 

m  -  g(hT}x)x  +  (2a  +  fi)hhxtixtt  +  ah^nxxtt  -  gh^(?ct  +  2  +  fi)hxrixxx 
-gh^^a  +  ^^Tixxxx  =  0  (31) 

Carrying  the  calculation  forward  to  the  point  of  obtaining  a  shoaling  relation  results  in  a  slight  mod¬ 
ification  to  the  expression  F  in  equation  (28) 

F  =  gk-\- (2a P)aP'kh  -  gh^(5a  +  2  +  (32) 

Kaihatu  and  Kirby  (1998)  determined  the  free  parameters  a  (for  dispersion)  and  p  (for  shoaling) 
using  a  least  squares  optimization  integrated  as  a  function  of  h/Lo.  Two  sets  of  parameters  were 
found.  The  first  set  optimized  the  shoaling  while  using  the  optimum  a  determined  by  Chen  and 
Liu  (1995).  This  set  (a  =  ->0.3855,  fi  =  -0.3540)  was  known  as  the  “dispersion  optimized”  (DO) 
set.  The  second  parameter  set  was  determined  by  finding  the  values  of  a  and  p  which  minimized 
the  global  error  in  (a,fi)  parameter  space.  This  set  (a  =  —0.4111,^  =  -0.3188)  was  denoted  the 
“dispersion  and  shoaling  optimized”  (DSO)  parameter  set.  It  is  noted  that  the  DSO  parameter  value 
of  O'  =  —0.4111  is  very  close  to  Witting’s  (1984)  optimum  dispersion  parameter  value  (found  via 
Padd  approximants)  of  a  =  —2/5. 


3.2.  Frequency  Domain  IVansformation  of  the  Equations  of 
Nwogu  (1993):  Nonlinear  Parabolic  Model 

If  nonlinearity  and  two-dimensionality  had  been  retained  when  deriving  equation  (31),  the  result 
would  have  been  (Kaihatu  and  Kirby,  1998) 

TJtt  -  •  (hVrj)  -  hV  •  [{Ua  •  V)Ua]  +  (2a  +  P)hVh  •  Vr]n  +  +  V  •  (??u«)r 

^  ghhSa  +2  +  p)Vh  .  ViV'^Ti)  -  =  0  (33) 

Complete  elimination  of  the  velocity  requires  assumption  of  time  periodicity  for  both  rj  and  u^. 
To  eliminate  Ua  from  the  nonlinear  terms,  we  can  use  the  time-periodic  form  of  equation  (18) 

(34) 

no) 

leading  to  the  time-periodic  equation  for  7)„.  To  facilitate  convenient  numerical  treatment,  we  make 
use  of  the  parabolic  approximation,  first  developed  by  Radder  (1979)  and  Lozano  and  Liu  (1980)  for 
the  linear  mild-slope  equation.  We  first  make  an  explicit  assumption  that  the  waves  are  propagating 
forward 

(35) 

The  complex  amplitudes  An(x,  y)  represent  phase-like  behavior  in  both  x  and  y.  The  phase  function 
is  integrated  only  in  x;  this  places  all  phase-like  behavior  in  y  in  the  complex  amplitude  An,  while 
allowing  explicit  slow  and  fast  variations  in  x.  The  consequence  is  that  the  angle  between  the  incident 
wave  and  the  x  direction  of  the  grid  remains  small  in  order  to  maintain  slowly-varying  wave-like 
behavior  of  An  in  the  y  direction. 

Because  of  the  third  and  fourth  derivatives  of  t]  present  in  equation  (33),  terms  proportional  to 
KxAnx,  Anxxy  and  other  higher-order  derivatives  will  be  generated.  We  thus  keep  a  higher  degree  of 
modulation  in  the  y  direction  than  in  the  x  direction  (following  the  ordering  of  Liu  et  al.,  1985),  lead¬ 
ing  to  a  parabolic  evolution  equation  for  An .  However,  because  the  phase  function  in  equation  (35) 
is  integrated  only  in  x,  but  k  is  a  function  of  both  jc  and  y ,  a  redefinition  is  required 

An  =ane^  ^ ix,y)dx 

where  kn  is  a  reference  wave  number  that  is  the  result  of  averaging  along  the  y  direction.  With  this 
redefinition,  the  model  equation  reads 

2i  ^hkn  +  rP'Qp’aknh^  -  j  anx 

-  2  ^hkn  +  rP’op-akniP'  -  Igh^k^  -f-  j  (fe„  -  k„)a„ 

+  i\.gkn  +  n^oP’ilu  -f  P)knh  -  gh^k^{5a  +  2-1-  P)'\hxan 
+  i  +  n^oP'ah^  -  6gh^kl  j  knxOn 

+  [5  +  n^o?-{2a  +  P)h  -  gh^k^iSa  +  2  +  /3)j  hyOny 
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n-\  ^  ^  ^  N-n  _  \ 

+2  J2  5afa^+/^^’/(^«+/-^/-^n¥A:  ^q  (37) 

,/=!  /=1  7 


where 

/(«-/)  +  gh 

(38) 

o  r#  I  ,  .2-, 

/(«  +  0a>2^  '  '  W/)]+ 

(39) 

This  is  the  primary  result  from  Kaihatu  and  Kirby  (1998). 

Kaihatu  (1994)  noted  that  the  ambiguity  which  affected  the  substitution  process  in  the  linear  terms 
also  appears  in  the  formulation  of  the  nonlinear  terms,  to  the  effect  that  several  different  nonlinear 
coefficient  sets  are  possible.  Unlike  the  linear  terms,  however,  there  are  no  corresponding  analytical 
metrics  for  determining  which  set  of  coefficients  have  the  most  desirable  properties,  outside  of  com¬ 
parisons  to  nonlinear  permanent  form  solutions.  The  general  complexity  of  the  Boussinesq  equations 
of  Nwogu  (1993),  particularly  in  retaining  shoaling  and  dispersive  properties  during  transformation 
into  the  frequency  domain,  has  motivated  investigations  into  developing  simpler  forms  of  the  equa¬ 
tions.  This  is  explored  in  the  next  section. 


3.3.  Frequency  Domain  TVansformation  of  the  Equations  of  Chen  and  Liu  (1995) 

Chen  and  Liu  (1995)  and  Kaihatu  and  Kirby  (1994)  investigated  using  the  extended  Boussinesq 
equations  of  Nwogu  (1993)  in  the  form  of  velocity  potential  (rather  than  velocity)  and  free  surface 
elevation.  The  resulting  equations  would  be  an  analogue  to  the  Boussinesq  equations  of  Wu  (1981), 
which  were  a  (<^,  t?)  form  of  the  classical  Boussinesq  equations  of  Peregrine  (1967).  The  use  of  ve¬ 
locity  potential  as  a  dependent  variable  simplifies  the  treatment  of  the  extended  Boussinesq  equations 
considerably,  since  ^  is  a  scalar  quantity. 

Chen  and  Liu  (1995)  began  with  the  boundary  value  problem  for  water  waves,  with  surface  bound¬ 
ary  conditions  scaled  for  shallow  water  waves  and  expanded  to  0(SypL^),  The  derivation  is  similar 
to  that  seen  in  Mei  (1983),  except  that  the  velocity  potential  is  taken  at  an  unknown  level  in  the 
water  column  rather  than  averaged  over  depth.  Kaihatu  and  Kirby  (1994)  expanded  the  equations  to 
0(8yfi^y  5/x^),  thereby  including  dispersive  effects  in  the  nonlinear  terms.  The  resulting  equations 
are 


»?»  +  V  .  [(/t  +  „)  +  V  •  V  jz„  V .  +  y  j 


+y  V  [V  •  (A V^a)]  -  y  vv2^„  I  +  T  =  0 


1  T 

+  2  (V^a)  -I- 


ZaV  • 


^at)  +  J 


+  A=0 


(40) 


(41) 


where  (pa  is  the  velocity  potential  at  Za-  The  terms  T  and  A  represent  terms  of  D(5/x^),  which  are 
zero  in  Chen  and  Liu  (1995).  In  Kaihatu  and  Kirby  (1994)  they  are 


T  =  V  .  (ijV^a)  -h  V  .  (|?V  V  .  (hV<Pa)]) 


+  TiV 
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A  =  -t,V-  {hV4>at)  +  V^„  .  [V2„  V  •  (hV<l>a)]  +  Za  (V^a  •  V)  (V  •  (A V^„)) 

+  •  (za  Vz„  V2^„)  +  ^  V^«z^  •  V(V2^a)  +  iv  ■  (/.V^„)  (43) 

Temporal  periodicity  was  assumed  for  the  velocity  potential  and  the  following  equation  developed 
for  the  amplitudes  of  the  velocity  potential  bn  (Chen  and  Liu,  1995,  although  expressed  somewhat 
differently) 


^  -  2  («  + 1)  1,3,2  j  +  21  (/.  +  -  2fc>3  +  1  y 

j^l  +  -  (1  +5a  +  vT+2a)Jfc2/,2j 

+  i  +  5a  +  VH^)]  j 

-  2  ^k„  -  2*3,3  ^  1  ^  j 

,  A-1  ^  _  _  N-n  _  _  _  \ 

—  —  I  "^Rbibn-ie^ ) 

^  \/=l  /=1  / 


where 


R  ==  4-  2nkikn-i  +  (n  -  l)kf  -  ah^  (^Ikfkn-i  +  nkf  +  (n  - 


(44) 


(45) 


S  =  {n+  l)kf  -  2nkik„+i  -  Ikj+i  +  ah^  ((«  +  /)*/*3^,  -  nkfkl^,  -  lk„+ikf)  (46) 

From  the  second-order  dynamic  free  surface  boundary  condition,  the  nonlinear  relation  between  the 
amplitudes  of  <pa  and  those  of  the  free  surface  rj  can  be  derived 


an 


i.2t  ,  '®n2at  t  2(0„0lh\n  ^  (OnZakn  ,  , 

— - an  Onyy  H - Hyt>ny - Onjc - njcOn 

8  8  8  8 


^^knxbn  +  —[1  -oik'^h^-2ahh„(k„-kn)]bn 
8  8 


1 


-  _  -  N-n 

f{ki+k„-i-k„)dx  _  2  ^  s'bib„+ie'S(ktH:l 
1=1 


(47) 
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where 

7 

II 

(48) 

s' = 

(49) 

The  extension  to  0(8fi^)  results  in  the  nonlinear  coefficients  (Kaihatu  and  Kirby,  1994) 


R  =  R-oih^  ((n  - l)kf  +{2n- l)kfk„-i  +  («  + 

_  „a2  ( (n-l)^+l(n-l)  +  fi\  2  2 
[  Kn-l)  p^n-l 


5  =  S  +  a/,2  {lk*^^  +  (2n+  l)kfk„+i  +  («  -  -  («  +  l)kf) 

+  nh^  {  +  0  +  .2^2 

S.s'-ahk  »? 


2jfe2 


(50) 


(51) 


(52) 


(53) 


The  0(5/x^)  terms  generally  tend  to  induce  a  slower  growth  in  the  generated  higher  harmonics 
in  shoaling  wave  applications.  Truncation  of  the  extended  Boussinesq  model  to  (with  all 

nonlinearity  retained)  was  performed  by  Wei  et  al.  (1995).  We  note  that  the  ambiguities  which  arose 
in  the  previous  frequency  domain  treatment  of  the  equations  of  Nwogu  (1993)  are  not  present  here. 


3*4*  Model  Evaluation 

We  noted  earlier  that  a  primary  motivation  for  development  of  the  extended  Boussinesq  models 
was  to  impart  linear  properties  which  mimicked  those  of  fiiUy-dispersive  linear  theory  for  a  wide 
range  of  water  depths,  thus  removing  a  substantial  obstacle  to  general  model  application.  In  this 
section  we  examine  how  well  the  resulting  frequency  domain  models  capture  wave  shoaling,  having 
insured  (via  the  substitution  process)  that  linear  dispersion  is  well  represented.  Additionally,  because 
the  resulting  models  are  parabolic,  we  will  also  examine  the  wide-angle  behavior  of  the  equations. 


3.4.  L  Wave  Shoaling 

As  mentioned  previously,  one  of  the  consequences  of  frequency  domain  transformation  of  the 
equations  of  Nwogu  (1993)  has  been  the  lack  of  any  guarantee  that  the  characteristics  of  the  original 
equations  survive  the  transformation  process.  Great  care  had  to  be  exercised  in  the  combination  of 
the  equations  such  that  the  advantageous  linear  dispersion  properties  could  be  maintained.  Additional 
steps  were  required  for  handling  shoaling.  In  contrast,  the  equations  of  Chen  and  Liu  (1995)  were 
relatively  simple  to  transform  into  the  frequency  domain,  particularly  with  respect  to  retaining  the 
dispersion  relation  of  the  equations  of  Nwogu  (1993). 
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Madsen  and  S0rensen  (1992)  noted  that  the  most  reliable  measure  of  the  effectiveness  of  the 
shoaling  mechanism  was  the  shoaling  gradient 


\Mx  ^  hx 
\A\  h  ^ 


(54) 


where  y  is  the  shoaling  gradient.  While  deviation  of  the  shoaling  gradient  from  linear  theoiy  may 
exaggerate  that  seen  in  the  shoaling  model  itself,  particularly  in  intermediate  water  depth  (Chen  and 
Liu,  1995),  it  is  a  convenient  measure  since  integration  is  not  required.  The  resulting  expressions 
for  y  for  the  shoaling  models  compared  herein  can  be  found  in  the  original  publications.  Kaihatu 
and  Kirby  (1998),  as  mentioned  previously,  used  two  sets  of  parameters:  the  DO  parameters  (a  = 
—0.3855,  =  —0.3540)  and  the  DSO  parameters  (cc  =  —0.4111,  p  =  —0.3188).  To  demonstrate  the 
effect  of  the  P  term,  we  also  include  the  (a  =  —0.3855,  ^  =  0)  case;  this  is  what  would  result  if  the 
ambiguity  in  the  substitution  process  detailed  earlier  had  not  been  realized.  The  shoaling  mechanisms 
of  Madsen  and  S0rensen  (1992)  and  Chen  and  Liu  (1995)  will  be  used  with  the  optimized  free 
parameters  determined  by  the  authors  in  the  original  publications.  We  note  that  the  model  of  Kaihatu 
and  Kirby  (1994)  shares  the  same  shoaling  characteristics  as  that  of  Chen  and  Liu  (1995). 

As  a  benchmark,  the  shoaling  gradient  from  fully-dispersive  linear  theory  is  used  (Madsen  and 
S0rensen,  1992) 


G'(l  +  |G'(l-cosh2it/j)) 

- (UW - 

where 
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sirihlkh 
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Figure  1  shows  a  comparison  between  the  different  shoaling  mechanisms  and  that  of  linear  theory.  Of 
the  five  shoaling  mechanisms,  those  of  Madsen  and  S0rensen  (1992)  and  Kaihatu  and  Kirby  (1998) 
with  DSO  parameters  compare  the  best,  with  a  slight  improvement  yielded  by  the  free  parameteriza¬ 
tion  of  the  latter  model.  In  contrast,  the  shoaling  mechanisms  of  Chen  and  Liu  (1995)  and  Kaihatu 
and  Kirby  (1998)  with  ^  =  0  compare  poorly,  particularly  the  latter  model.  The  model  of  Kaihatu 
and  Kirby  (1998)  may  represent  the  limit  of  optimum  shoaling  performance  possible  with  two  free 
parameters.  Sch^er  and  Madsen  (1995),  using  a  Pad6  [4, 4]  approximation  and  more  free  parame¬ 
ters,  developed  an  expression  for  y  which  exhibited  virtually  no  deviation  from  linear  theory  for  the 
complete  range  of  water  depths. 


5.4.2.  Wide-Angle  Behavior  of  the  Parabolic  Equations 

One  consequence  of  the  parabolic  approximation  is  the  assumption  that  the  angle  between  the 
x-coordinate  of  the  grid  and  the  incident  wave  direction  remains  small.  Methods  for  ameliorating 
this  problem  have  generally  taken  the  form  of  retention  of  higher-order  derivative  terms  with  coeffi¬ 
cients  which  can  be  determined  by  Pad6  approximations  (Booij,  1981;  Kirby,  1986a)  or  by  rational 
approximations  (Kirby,  1986b),  and  have  been  shown  to  work  well  for  parabolic  approximations  of 
the  raild-slope  equation  (Berkhoff,  1972;  Smith  and  Sprinks,  1975).  The  suitability  for  parabolic  ap¬ 
proximations  of  the  Boussinesq  equations,  however,  have  not  been  established  except  by  comparisons 
with  data  (for  example,  Kirby,  1990).  In  this  section  we  analyze  the  effectiveness  of  the  parabolic 
approximation  used  in  the  development  of  the  frequency  domain  extended  Boussinesq  models.  We 
note  that  the  models  of  Chen  and  Liu  (1995)  and  Kaihatu  and  Kirby  (1994;  1998)  reduce  to  the  same 
basic  parabolic  form  in  the  linear  limit. 
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Numerical  treatment  of  parabolic  wave  models  generally  involves  the  Crank-Nicholson  discretiza¬ 
tion  scheme,  which  is  second-order  accurate  in  both  horizontal  directions.  A  tridiagonal  matrix  is 
formed  at  each  x  location  and  is  usually  solved  using  a  Thomas  algorithm.  The  highest  order  deriv¬ 
ative  of  A  which  would  allow  this  treatment  is  Axyy  (Kirby,  1986a),  which  was  not  retained  in  our 
previous  development.  Expanding  the  time-periodic  form  of  equation  (33)  into  its  horizontal  com¬ 
ponents,  neglecting  nonlinear  and  bottom  slope  terms,  substituting  in  equation  (35)  and  retaining  the 
Axyy  term  that  results  yields 

2j  +  Q^kh^  -  2g{khf  j  Ax  +  -  Igh^k^  ^ j  Ayy 

+  +  (57) 

We  note  that  many  terms  were  truncated  in  the  substitution  process  used  to  derive  equation  (57), 
specifically  third  and  fourth  derivatives  with  respect  to  y.  Additionally,  similar  derivatives  with  re¬ 
spect  to  X  are  represented  only  in  the  differentiation  of  the  phase  function,  generating  terms  propor¬ 
tional  to  ,  To  ascertain  the  effect  this  truncation  has  on  the  accuracy  of  wave  propagation,  we  need 
to  transform  the  equation  back  into  ?}  by 

A  =  (58) 

leading  to 

2i  -f  oP’kh^a  -  Igk^h^  -f  j  r}x  -h  ct?kh^ot  ~  Igk^h^  -h  j  rj 
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+  \^h  +  oP'h^a  -  Igh^k^  O]  '*‘^*^*  (“  +  j) 

+  4k‘^gh^(a  +  ^'^fiyy=0  (59) 

Making  a  final  substitution 

$  =  (60) 


where  k;,;  and  ky  denote  wave  number  vector  components  in  the  x  and  y  directions,  respectively,  we 
can  obtain  an  expression  for  ky  in  terms  of  k  and  kx 


ky  = 


2{k  —  kx)  -f  Qp'kh'^ct  —  Igk^h^  1 

gh 

AhHk'^-kkx)\ 

1  + 

(61) 


The  neglect  of  the  Axyy  term  would  have  the  effect  of  zeroing  out  the  kkx  term  in  the  denominator 
of  equation  (61).  The  baseline  for  comparison  is  the  circle 


(62) 


which  was  derived  by  substitution  of  equation  (60)  into  the  Helmholtz  equation.  Additionally,  we 
also  compare  to  the  small-angle  parabolic  approximations  used  in  the  mild-slope  equation  (Kirby, 
1986b) 


(63) 


Figure  2.  Analysis  of  the  parabolic  approximations  used  for  various  models.  Relation  from  elliptic  model  is 
benchmark. 
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and  the  wide  angle  expression  from  the  Fade  approximation  of  Kirby  (1986a) 


4*2 

J 

(64) 


Figure  2  shows  a  comparison  of  the  above  expressions  relating  kx  to  ky  At  appears  that  the  parabolic 
models  of  the  extended  Boussinesq  equation  have  slightly  worse  characteristics  at  oblique  angles 
than  the  small-angle  approximation  of  the  mild-slope  equation,  with  the  Axyy  term  imparting  no 
appreciable  improvement.  This  could  be  due  to  the  considerable  amount  of  information  contained 
in  the  terms  (among  others)  discarded  when  the  form  (equation  (35)  is  substituted  and  the 

parabolic  approximation  made.  Optimizations  similar  to  Kirby’s  (1986a;  1986b)  developments  with 
the  mild-slope  equation  could  be  implemented  here. 


4.  NONLINEAR  MILD-SLOPE  EQUATION  MODELS 

An  alternative  approach  to  developing  nonlinear  frequency  domain  models  involves  incorporating 
nonlinear  effects  into  models  already  equipped  with  fully-dispersive  transformation  characteristics. 
The  mild-slope  equation  (Berkhoff,  1972;  Smith  and  Sprinks,  1975)  simulates  wave  refraction,  shoal¬ 
ing  and  diffraction  over  mildly- varying  bathymetry;  its  applicability  has  been  greatly  increased  with 
the  advent  of  the  parabolic  approximation  (Radder,  1979;  Lozano  and  Liu,  1980)  encountered  in 
earlier  sections. 

Bryant  (1973,  1974)  first  studied  the  efficacy  of  developing  fully-dispersive  models  with  second- 
order  nonlinear  characteristics.  He  developed  a  spatially-periodic  solution  to  the  truncated  Laplace 
boundary  value  problem  (thereby  retaining  full  frequency  dispersion)  and  compared  numerical  eval¬ 
uations  of  this  solution  to  those  from  various  forms  of  the  Korteweg-de Vries  (KdV)equation.  Fur¬ 
thermore,  Bryant  (1974)  also  demonstrated  that  three  harmonic  amplitudes  of  his  spatially-periodic 
solution  matched  those  of  third-order  Stokes  waves,  as  did  the  nonlinear  dispersion  characteristics. 

Keller  (1988)  derived  coupled  nonlinear  equations  derived  from  the  shallow  water  equations,  the 
Boussinesq  equations  and  the  Euler  equations,  all  truncated  to  two  harmonics.  He  showed  that,  in  the 
shallow  water  limit,  the  equations  reduced  to  identical  forms. 

Agnon  et  al.  (1993)  developed  a  unidirectional  shoaling  model  based  on  a  nonlinear  extension 
of  fully-dispersive  linear  shoaling.  The  linear  part  of  the  resulting  model  contained  fully-dispersive 
shoaling,  and  triad  interactions  between  wave  components  described  the  nonlinear  evolution.  This 
was  later  extended  to  two-dimensional  propagation  by  Kaihatu  and  Kirby  (1995),  Tang  and  Ouelette 
(1997)  (both  parabolic  models)  and  Agnon  and  Sheremet  (1997)  (hyperbolic  model).  The  paper  by 
Kaihatu  and  Kirby  (1995)  also  detailed  the  inclusion  of  a  surf  zone  dissipation  mechanism.  Tang 
and  Ouelette  (1997)  extended  the  model  of  Kaihatu  and  Kirby  (1995)  by  including  diffraction  and 
bottom  slope  effects  in  the  nonlinear  terms  of  their  model. 

To  explain  some  of  the  subtle  points  outlined  later  in  this  section,  we  briefly  outline  the  develop¬ 
ment  of  the  nonlinear  mild-slope  equation  described  by  Kaihatu  and  Kirby  (1995).  We  start  from  the 
boundary  value  problem  for  water  waves,  with  free  surface  conditions  expanded  to  second-order  in 
wave  amplitude.  We  first  assume  that  the  solution  can  be  expressed  as  a  superposition 

N 

Hx,  y,  Z,  f)  =  fn(z)$ni^>  y)e~“^*  +  C.C. 
n=l 


(65) 
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where  is  complex,  and 
-  ^  cosh^n(^  +  g) 

"  coshknh 

where  the  dispersion  relation  is  that  of  linear  theory 
co^  =  gkn  tanh^n/i 

We  then  make  use  of  Green’s  Second  Identity  on  the  variables  /„  and  and  then  use  resonant  triad 
interaction  theory  to  create  a  time-periodic  evolution  equation  for  We  then  assume  a  propagating 
wave 

0„=-— (68) 

(On 

We  substitute  equation  (68)  and  its  conjugate  into  the  time-periodic  equation  for  and  employ  the 
parabolic  approximation  to  justify  neglecting  d^Anldx^  terms  in  the  resulting  equation.  We  then 
make  use  of  a  phase  function  redefinition  similar  to  that  done  for  the  extended  Boussinesq  models  in 
earlier  sections.  This  results  in 
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where 
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(69) 
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z  =  — ^ — [a>lkikn+i  +  {kn+i  -  */)(&)„+/*/  +  o)ik„+i)(o„] 
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-  y  io>f  -  (Oi(OnM  +  (of^^{)  (71) 

Equation  (69)  comprises  the  model  of  Kaihatu  and  Kirby  (1995). 

Equation  (68)  is  derived  from  the  first-order  dynamic  free  surface  boundary  condition 

<t>f^gri  =  0\  z-0  (72) 

and  is  a  transformation  from  amplitudes  of  velocity  potential  to  those  of  the  free-surface  elevation. 
This  transformation  is  first-order,  and  does  not  include  the  nonlinear  terms  inherent  in  the  dynamic 
free  surface  boundary  condition.  Eldeberky  and  Madsen  (1999)  determined  that  the  neglect  of  these 
second-order  terms  had  the  effect  of  underpredicting  the  nonlinear  energy  transfer,  particularly  with 
respect  to  the  superharmonic  energy  transfer.  They  used  successive  approximations  to  invert  the 
second-order  dynamic  free  surface  boundary  condition  and  eliminate  the  velocity  potential,  resulting 
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in  an  evolution  equation  that  could  be  conveniently  solved  for  the  free  surface  amplitudes  alone.  The 
linear  terms  are  the  same  as  those  of  earlier  models,  but  the  nonlinear  coefficients  are 

Y  =  Y  Q)nkikfi—i  H - —(o)f  -jrcoicofi—i  +co^_i)  (73) 

coicOn^l  g 

+  w^+;)  (74) 

where 

p+  _  -kn)Cgn 

(On  ^  ^ 


_  2{kn^l  —ki—  kn)Cgn 
- 

Eldeberky  and  Madsen  (1999)  demonstrated  that  the  above  terms  contributed  substantially  to  the 
superharraonic  energy  transfer.  We  note  here  that  Tang  and  Ouelette  (1997)  also  took  the  second- 
order  terms  in  this  transformation  into  account,  though  in  a  manner  different  than  Eldeberky  and 
Madsen  (1999). 

Kaihatu  (2001)  took  a  different  approach  and  derived  the  required  correction  to  make  the  result¬ 
ing  equations  consistent.  This  correction  was  derived  from  the  second-order  dynamic  free  surface 
boundaiy  condition,  and  is  applied  whenever  the  free  surface  is  required 

1  f /•  “  -  -  \ 

a„=a„  +  —  I  Y"aia„-ie'  f ik,+k„.,-kn)dx  +  2j2  Z"afa„+i^  f(k„+,-ki-k„)dx  j 

where  an  is  the  total  amplitude  of  the  free  surface  to  second-order,  %  is  the  solution  to  equation  (69), 
and 


Y"  =  co}  +  coia>„_,  +  o>l_i  -  8^^^  (78) 

Z"  =  o)f-  coi<o„+i  +  .^^±L 

Kaihatu  (2001)  also  investigated  the  effect  of  this  correction  by  comparing  permanent  form  solutions 
of  the  one-dimensional  version  of  equation  (69),  with  and  without  the  correction  (equation  (77)),  to 
third-order  Stokes  theory.  While  the  correction  did  not  noticeably  enhance  phase  speed  comparisons 
to  the  theory,  it  did  improve  comparisons  of  the  respective  free  surfaces.  Additionally,  Kaihatu  (2001) 
also  extended  the  parabolic  model  to  include  wide-angle  propagation  terms,  which  notably  improved 
model  performance  when  compared  to  laboratory  measurements  of  waves  propagating  over  a  shoal 
(Chawla  et  al,  1998).  These  wide-angle  propagation  terms  are  essentially  those  of  Kirby  (1986a)  and 
therefore  share  the  same  accuracy  in  oblique-angle  propagation  as  equation  (64). 
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5.  BREAKING  WAVE  MODELS 


For  general  utility  in  solving  nearshore  wave  propagation  problems,  consideration  of  wave  break¬ 
ing  and  surf  zone  decay  is  required.  Due  to  the  nature  of  the  evolution  equations,  such  breaking  and 
decay  descriptions  must  necessarily  be  statistical.  Two  formulations  are  generally  used  for  this  ap¬ 
plication;  those  of  Battjes  and  Janssen  (1978)  and  Thornton  and  Guza  (1983),  though  several  others 
are  extant  (for  example.  Dally,  1990;  Battjes  and  Groenendijk,  2000). 

Battjes  and  Janssen  (1978)  assumed  that  the  probability  distribution  of  breaking  waves  could  be 
described  as  a  Rayleigh  distribution,  where  the  percentage  of  breaking  waves  at  a  particular  location 
is  related  to  the  area  under  the  truncated  probability  distribution.  This  percentage  Qi,  was  determined 
by  solution  of  the  following  implicit  relation 

l-gfc 

Ine*  \Hmax) 

where  Hrms  is  the  root-mean-square  wave  height  and  Hmax  the  maximum  wave  height  in  the  distri¬ 
bution.  Framing  the  energy  dissipation  from  breaking  waves  as  an  energy  balance 

=  -eb  (81) 


where  E  is  wave  energy,  Battjes  and  Janssen  (1978)  determined  an  expression  for  the  energy  dissi¬ 
pation 

1  -  o 

—  '^Qbf  P  syntax  (82) 

where  /  is  the  average  frequency  of  the  spectrum.  The  maximum  wave  height  Umax  is  determined 
by 


(83) 


where  k  is  the  average  wave  number.  An  expression  for  y  based  on  the  deep  water  wave  steepness 
was  found  by  Battjes  and  Stive  (1985). 

Thornton  and  Guza  (1983)  extended  the  Battjes  and  Janssen  (1978)  model  by  accounting  for  the 
transformation  of  the  wave  height  probability  distribution  through  the  surf  zone.  They  hypothesized 
a  distribution  of  breaking  waves  as  being  a  weighted  Rayleigh  distribution  with  tunable  parameters. 
The  energy  dissipation  for  a  single  bore  was  then  integrated  through  this  probability  distribution  to 
obtain 


_  3^ 

*  7^*5 


(84) 


where  B  and  y  are  free  parameters.  Reasonable  fit  to  field  data  was  found  with  y  =  0.42  and 
B  =  1.3  —  1.7  (Thornton  and  Guza,  1983).  Mase  and  Kirby  (1992)  incorporated  this  dissipation 
mechanism  into  their  “hybrid  KdV”  shoaling  model,  which  used  full  linear  shoaling  with  shallow 
water  nonlinearity.  This  particular  study  also  detailed  a  laboratory  experiment  in  which  a  spectrum 
of  waves  was  allowed  to  shoal  and  break  over  a  long  sloping  bottom.  The  unique  feature  of  this 
experiment  is  that  the  energy  at  the  peak  of  the  spectrum  is  in  intermediate  water  depth  at  the  wave 
maker;  this  would  serve  as  a  severe  test  of  shoaling  models,  and  would  invalidate  those  limited  to 
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weak  dispersion  (for  example,  models  based  on  the  classical  Boussinesq  equations).  This  dissipation 
mechanism  was  also  included  in  the  linear  spectral  parabolic  model  of  Chawla  et  al.  (1998),  as  well 
as  the  nonlinear  shoaling  model  of  Kaihatu  and  Kirby  (1995). 

Both  Battjes  and  Janssen  (1978)  and  Thornton  and  Guza  (1983)  developed  their  dissipation  mech¬ 
anisms  purely  as  lumped  parameter  models,  dependent  only  on  integrated  properties  of  the  spectrum 
with  no  other  details.  For  use  in  phase-resolving  frequency  domain  models  such  as  those  detailed  in 
this  study,  some  assumptions  concerning  the  distribution  of  the  dissipation  over  the  frequency  range 
must  be  made.  The  majority  of  nonlinear  models  of  this  type  assume  an  equal  weighting  of  dissipa¬ 
tion  across  the  entire  frequency  range  (Eldeberky  and  Battjes,  1996;  Eldeberky  and  Madsen,  1999); 
this  assumption  is  also  used  in  linear  spectral  models  (Chawla  et  al.  1998).  Alternatively,  Mase  and 
Kirby  (1992),  Kaihatu  and  Kirby  (1995)  and  Kirby  and  Kaihatu  (1996)  used  a  distribution  which  as¬ 
sumes  a  frequency-squared  weighting  of  the  dissipation  term,  thus  accounting  for  nonlinear  transfer 
of  energy  from  low  to  high  frequencies  due  to  triad  interactions.  Chen  et  al.  (1997)  demonstrated 
(using  the  model  of  Chen  and  Liu  1995)  that,  while  the  frequency-squared  weighting  did  not  affect 
the  spectral  shape  significantly,  it  did  offer  greater  accuracy  in  estimating  skewness  and  asymmetry. 
Both  quantities  are  indications  of  wave  shape;  for  example,  negative  asymmetry  corresponds  to  a 
forward-pitched  shape  of  the  wave  field,  redolent  of  surf  zone  waves. 

The  dissipation  mechanisms  are  usually  formulated  in  terms  of  energy  E,  Some  manipulation  is 
required  to  implement  these  into  complex-amplitude  evolution  equations.  Mase  and  Kirby  (1992) 
implemented  the  Thornton  and  Guza  (1983)  dissipation  by  starting  with  the  conservation  of  energy 
flux  equation  with  a  damping  term  added 

A„^  +  ^^A„+oiA„=0  (85) 

where  a  is  the  dissipation  coefficient.  Multiplying  this  equation  by  the  conjugate  amplitude,  adding 
it  to  its  conjugate  equation,  then  summing  over  all  components,  we  obtain 

N  N 

J^{Cgn\An\^)^  =  -2j2^nCgn\An\^  (86) 
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Then,  assuming  shallow  water  and  switching  to  an  energy  definition,  we  obtain 


Pg^fgh 
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where  p  is  mass  density  of  water.  Equating  this  to  the  dissipation  function  in  equation  (84)  yields 
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We  now  require  a  frequency  distribution  for  otn-  Mase  and  Kirby  (1992)  investigated  the  trends 
evidenced  in  a  back-calculation  of  Un  from  the  data,  and  determined  that  a  strong  dependence  for 
the  dissipation  existed,  analogous  to  a  frequency  domain  transformation  of  the  Burgher’s  equation 
with  viscous  damping.  A  reasonable  representation  of  this  frequency  dependence  can  be  achieved  by 
assuming  the  following  form  for 


=  tt/iO  + 
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(89) 
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where 

a„o  =  Fp(x)  (90) 


«nl  = 


ELlfi\An\^ 


(91) 


This  essentially  splits  the  frequency  dependence  into  a  balance  between  one  that  is  flat  across  the 
frequency  range  (5^0)  that  is  weighted  to  the  square  of  the  frequency  This  split  is 

controlled  by  the  parameter  F.  Mase  and  Kirby  (1992)  found  that  F  =  0.5  seemed  to  work  best  for 
their  experimental  data,  though  Chen  et  al.  (1997)  noted  that  F  =  0  worked  best  if  the  comparison  at 
the  shallowest  gauge  of  the  experiment  {h  =  .025  m)  was  neglected. 

Eldeberky  and  Battjes  (1996)  and  Eldeberky  and  Madsen  (1999)  both  used  the  dissipation  of  Bat- 
tjes  and  Janssen  (1978),  and  assumed  that  the  dissipation  was  weighted  equally  across  all  frequencies. 
By  realizing  that  amplitude  changes  are  half  those  of  changes  in  energy,  they  found  that 
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where  the  denominator  represents  the  total  energy  flux  in  the  spectrum.  Eldeberky  and  Madsen  (1999) 
show  comparisons  of  their  model  (with  the  inclusion  of  the  second-order  relation  between  (p  and  rj 
discussed  earlier)  with  the  data  of  Mase  and  Kirby  (1992).  While  they  showed  improved  comparisons 
of  skewness  relative  to  the  original  model  of  Kaihatu  and  Kirby  (1995),  they  were  unable  to  repli¬ 
cate  the  increase  in  negative  asymmetry  seen  in  the  data.  Kaihatu  (2001),  using  the  dissipation 
weighting  detailed  in  Mase  and  Kirby  (1992)  as  well  as  the  second-order  transformation  between  0 
and  demonstrated  improved  asymmetry  predictions  relative  to  those  seen  in  Eldeberky  and  Madsen 
(1999)  by  using  F  =  0  in  equation  (90). 


6.  COMPARISONS  TO  DATA 
6.1.  Whalin(1971) 

In  this  section  we  compare  several  of  the  described  models  to  the  experimental  data  of  Whalin 
(1971),  who  conducted  a  laboratory  study  to  investigate  the  limits  of  linear  refractive  wave  propa¬ 
gation  theory.  He  generated  sinusoidal  waves  with  periods  of  1,  2  and  3  seconds  and  ran  them  over 
bathymetry  resembling  a  tilted  cylinder.  The  plan  view  of  the  experimental  layout  is  shown  in  Fig¬ 
ure  3.  Wave  gauges  located  down  the  centerline  of  the  tank  measured  the  free  surface  elevation;  these 
measurements  were  then  processed  to  yield  wave  harmonic  amplitudes.  The  experimental  conditions 
ranged  from  deep  (/i^  ^  2)  to  shallow  (fj?"  ^  0.2)  water  at  the  wave  maker.  We  will  concentrate  on 
the  1  second  period  case.  Comparisons  to  other  cases  in  the  experiment  are  detailed  in  the  original 
papers. 

For  this  T  =  1  s  case  (ao  —  0.0195  m)  N  =  2  harmonics  were  used.  This  is  in  concert  with  the 
work  detailed  in  the  original  studies;  the  inclusion  of  additional  harmonics  did  not  make  a  signif¬ 
icant  difference.  For  each  case,  all  wave  energy  was  placed  in  the  first  harmonic,  with  the  higher 
harmonics  initialized  with  zero  energy.  We  compare  the  following  models:  the  extended  Boussi- 
nesq  model  of  Kaihatu  and  Kirby  (1998)  with  DSO  parameters  (a  -  -0.4111,  fi  =  -0.3188);  the 
extended  Boussinesq  model  (expressed  in  terms  of  (p  and  t])  of  Kaihatu  and  Kirby  (1994);  and 
the  nonlinear  mild-slope  equation  model  of  Kaihatu  and  Kirby  (1995)  with  the  second-order  cor¬ 
rection  of  Kaihatu  (2001).  The  exclusion  of  the  second-order  correction  of  Kaihatu  (2001)  did 
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not  yield  significant  difference  relative  to  inclusion.  Additionally,  comparisons  between  the  DO 
(a  ==  —0.3855,  p  =  —0.3540)  and  DSO  parameters  for  this  experiment  are  shown  in  the  original 
paper  of  Kaihatu  and  Kirby  (1998)  and  are  not  included  here.  Lastly,  the  model  of  Chen  and  Liu 
(1995)  was  not  used  in  these  comparisons;  they  used  a  second-order  bound  wave  solution  to  force 
their  model  for  the  7  =  1  s  case,  citing  a  large  phase  mismatch  between  the  free  and  bound  second 
harmonic  otherwise.  However,  we  are  more  interested  in  testing  the  model  in  general  wave  propaga¬ 
tion  scenarios,  where  a  priori  consideration  of  the  free  or  bound  wave  nature  of  the  forcing  would 
not  be  desirable.  Linear  wave  forcing  of  the  model  of  Chen  and  Liu  (1995)  for  this  case  reveals 
significant  oscillation  of  the  second  harmonic  over  the  domain. 

Comparisons  of  the  model  to  the  7  =  1  s,  =  0.0195  m  case  are  shown  in  Figure  4.  In  this 
case  the  model  of  Kaihatu  and  Kirby  (1995)  with  the  second-order  correction  (Kaihatu,  2001)  seems 
to  work  best  for  capturing  evolution  characteristics  of  both  harmonics.  The  extended  Boussinesq 
models  of  Kaihatu  and  Kirby  (1994;  1998)  appear  to  underpredict  the  energy  transfer  to  the  second 
harmonic  in  the  focal  region.  The  oscillations  in  the  second  harmonic  as  predicted  by  Kaihatu  and 
Kirby  (1994)  are  of  the  same  order  as  those  shown  in  Chen  and  Liu  (1995),  and  are  less  severe  than 
would  be  seen  if  the  Chen  and  Liu  (1995)  model  had  been  forced  with  a  linear  wave.  This  indicates 
that  the  presence  of  terms  appear  to  have  a  stabilizing  effect  on  the  model,  as  commented 

by  Tang  and  Ouelette  (1997).  We  also  note  here  that  the  model  of  Chen  and  Liu  (1995)  with  bound 
wave  forcing  agrees  with  the  second  harmonic  data  for  this  case  better  than  the  model  of  Kaihatu  and 
Kirby  (1994)  with  linear  wave  forcing.  Lastly,  Tang  and  Ouelette  (1997)  show  a  better  match  to  the 
7  =  1  s  case  than  seen  here,  possibly  due  to  the  inclusion  of  nonlinear  diffraction  terms. 

6.2.  Mase  and  Kirby  (1992) 

Mase  and  Kirby  (1992)  conducted  a  series  of  laboratory  experiments  in  which  irregular  waves 
(Pierson-Moskowitz  spectrum)  were  generated  and  allowed  to  shoal  and  break  over  a  slope.  The 
experimental  set-up  is  shown  in  Figure  5.  The  case  studied  here  had  a  peak  period  7p  =  1  s,  leading 
to  kh^  2  at  the  wave  maker.  Significant  wave  breaking  occurred  in  this  experiment  beginning  near 
the  gauge  zXh  —  Q.  175  m.  Time  series  of  free  surface  elevations  were  collected  at  20  Hz  and  divided 
into  7  realizations  of  2048  points  each.  Each  realization  was  then  put  into  a  Fast  Fourier  Transform 
(FFT);  the  resulting  energy  density  spectra  were  both  Bartlett-averaged  across  all  seven  realizations 
and  band-averaged  across  eight  neighboring  frequency  bands.  A  gauge  located  0.20  cm  offshore 
of  the  toe  of  the  slope  provided  the  initial  condition.  Tests  with  iV  =  300  and  yV  =  500  frequency 
components  were  run  through  the  models  for  each  realization,  leading  to  a  maximum  frequency  of 
3  Hz  and  5  Hz  respectively. 

We  compare  two  models  to  the  experimental  data:  the  nonlinear  mild-slope  equation  model  of  Kai¬ 
hatu  and  Kirby  (1995),  with  corrections  by  Kaihatu  (2001);  and  the  extended  Boussinesq  frequency 
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Figure  4.  Comparison  of  nonlinear  frequency  domain  models  data  of  Whalin  (1971).  Solid  line:  model  of  Kai- 
hatu  and  Kirby  (1998)  with  DSO  parameter!.  Dashed  line:  model  of  Kaihatu  and  Kirby  (1995)  with  corrections 
from  Kaihatu  (2001).  Dash>dot  line:  model  of  Kaihatu  and  Kirby  (1994).  Top:  first  harmonic  amplitudes.  Bottom: 
second  harmonic  amplitudes. 


domain  model  of  Kaihatu  and  Kirby  (1998)  with  the  DSO  parameters.  Both  models  were  equipped 
with  the  dissipation  mechanism  of  Thornton  and  Guza  (1983)  and  the  frequency-distribution  method¬ 
ology  detailed  by  Mase  and  Kirby  (1992).  Both  F  =  0.5  and  F  =  0  were  used,  the  latter  correspond¬ 
ing  to  a  full  frequency-squared  weighting  of  dissipation. 

Comparisons  of  spectra  at  various  locations  are  shown  in  Figure  6,  with  N  =  300  and  F  ==  0. 
It  is  apparent  that,  while  both  models  tend  to  compare  well  to  the  data  at  the  frequency  peak  and 
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Figure  6.  Comparisons  of  spectra  from  models  to  data  from  Mase  and  Kirby  (1992):  N  =  300  frequencies, 
F  =  0.  Solid  line:  data.  Dashed  line:  model  of  Kaihatu  and  Kirby  (1995)  with  correction  of  Kaihatu  (2001). 
Dash^dot  line:  model  of  Kaihatu  and  Kiiby  (1998)  with  DSO  parameters.  Top  left:  h  =  0.47  m.  Top  right: 
h  =  0.2  m.  Bottom  left:  h  =  0.125  m.  Bottom  right:  h  =  0.05  m. 


lower,  the  extended  Boussinesq  model  of  Kaihatu  and  Kirby  (1998)  significantly  underpredicts  the 
high  frequency  evolution  (/  >  1.75  Hz).  The  case  of  N  =  500  components  reveal  similar  trends.  No 
significant  differences  occur  between  the  F  =  0  and  F  =  0.5  cases. 
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Figure  7.  Comparisons  of  skewness  and  asymmetry  from  models  to  data  from  Mase  and  Kirby  (1992):  N  =  5Q0 
frequencies.  Open  circles:  data.  Solid  and  dashed  lines:  model  of  Kaihatu  and  Kirby  (1995)  with  correction  of 
Kaihatu  (2(K)1)  with  F  =  0.5  and  F  —  0,  respectively.  Dash-dot  and  dash-x  lines:  model  of  Kaihatu  and  Kirby 
(1998)  with  DSO  parameters  and  F  =  0.5  and  F  =:  0,  respectively.  Top:  skewness.  Bottom:  Negative  asymmetry. 


Differences  between  the  F  =  0  and  F  =  0.5  cases  begin  to  appear  when  comparing  skewness  and 
asymmetry.  Figure  7  shows  a  comparison  of  both  models,  each  using  F  =  0,5  and  F  =  0,  to  the  skew¬ 
ness  and  asymmetry  values  from  the  data  for  the  case  of  N  -  500.  (It  was  shown  by  Kaihatu,  2001 
that  good  comparisons  to  higher-order  moments  required  =  500  frequency  components.)  Both 
skewness  and  asymmetry  are  better  predicted  by  the  model  of  Kaihatu  and  Kirby  (1995),  with  cor¬ 
rections  by  Kmhatu  (2001),  than  with  the  extended  Boussinesq  model  of  Kaihatu  and  Kirby  (1998). 
Additionally,  F  =  0.5  results  in  a  slightly  better  skewness  prediction  and  a  slightly  worse  asymmetry 
prediction  than  F  =  0,  except  for  the  gauge  nearest  to  shore.  At  that  location  the  modeled  asymme¬ 
try  falls  off  dramatically  if  F  =  0  in  either  shoaling  model.  The  reason  for  this  sudden  dropoff  is 
not  clear.  We  note  here  that  this  behavior  at  the  last  gauge  is  not  an  indictment  of  the  F  =  0  value; 
Kennedy  et  al.  (2000)  showed  excellent  agreement  between  the  time  domain  Boussinesq  model  of 
Wei  et  al.  (1995)  and  the  Mase  and  Kirby  (1992)  data.  An  eddy  viscosity  mechanism  was  used  for 
dissipation  in  the  time  domain  model,  the  formulation  of  which  is  equivalent  to  F  =  0  in  the  fre¬ 
quency  domain.  Rather,  the  problem  in  the  frequency  domain  may  lie  in  the  use  of  the  bulk  energy 
dissipation  mechanism  used  here;  it  is  quite  likely  that  a  frequency  domain  formulation  of  the  local 
eddy  viscosity  mechanism  in  the  Kennedy  et  al.  (2000)  model  would  address  the  problems  seen  at 
the  shallowest  gauge. 
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6.3.  Discussion 

Based  on  the  data-model  comparisons  shown,  it  appears  that  frequency  domain  formulations  of 
the  extended  Boussinesq  model  of  Nwogu  (1993)  may  be  inferior  to  those  based  on  the  nonlinear 
mild-slope  equation  (Kaihatu  and  Kirby,  1995;  Kaihatu,  2001).  One  possible  reason  for  the  relatively 
poor  performance  of  the  frequency  domain  model  of  Kaihatu  and  Kirby  (1998)  may  lie  in  the  order 
of  truncation  of  the  original  equations  of  Nwogu  (1993).  These  equations  were  truncated  to  retain 
only  0(3,  ^^)  terms.  As  seen  in  the  (^,  ?/)  model  of  Kaihatu  and  Kirby  (1994),  it  is  possible  to  retain 
OifiiJ?’)  terms  and  still  treat  the  equations  using  the  frequency  domain  transformation;  perhaps  the 
inclusion  of  these  terms  in  the  frequency  domain  model  of  Kaihatu  and  Kirby  (1998)  would  improve 
the  results  (though,  in  the  author’s  experience,  the  derivation  process  is  exceedingly  complicated  and 
ambiguous). 

Interestingly,  however,  Wei  and  Kirby  (1995)  show  very  good  agreement  between  their  numer¬ 
ical  treatment  of  the  model  of  Nwogu  (1993)  and  the  Mase  and  Kirby  (1992)  data.  As  mentioned 
early  in  the  chapter,  there  is  often  some  discrepancy  between  the  performance  of  the  original  time 
domain  model  and  its  corresponding  frequency  domain  formulation.  Additionally,  Kofoed-Hansen 
and  Rasmussen  (1998)  showed  that  the  model  of  Madsen  and  Sdrensen  (1993)  (a  frequency  domain 
shoaling  model  based  on  the  extended  Boussinesq  model  of  Madsen  et  al.,  1991)  demonstrates  fa¬ 
vorable  agreement  with  the  Mase  and  Kirby  (1992)  data  set,  despite  the  order  of  truncation  being  the 
same  as  that  in  Kaihatu  and  Kirby  (1998),  This  is  likely  due  to  differences  in  the  underlying  time 
domain  equations  and  in  the  derivation  procedure  of  the  corresponding  frequency  domain  models. 


7.  STOCHASTIC  MODELS 

While  frequency  domain  models  offer  great  utility  in  simulating  shallow  water  processes,  this  is 
done  at  some  computational  expense,  particularly  when  they  are  initialized  with  smoothed  spectra 
(as  would  be  the  case  for  field  studies,  for  example).  Initial  phases  are  required  for  these  models, 
so  multiple  temporal  realizations  are  run  with  random  initial  phases  and  the  results  averaged  until 
acceptably  smooth,  which  can  require  considerable  computation  time.  This  has  motivated  research 
into  stochastic  models  of  triad  interactions.  In  this  section,  we  briefly  describe  a  few  approaches, 
referring  interested  readers  to  the  papers  cited  herein  and  the  overview  of  Agnon  and  Sheremet 
(2000)  for  further  details. 

Abreau  et  al.  (1992)  developed  a  statistical  model  for  triad  wave  evolution,  suitable  as  a  source 
function  in  a  spectral  balance  model  such  as  WAM  (Komen  et  al.,  1994)  or  SWAN  (Booij  et  al.,  1999). 
The  model  was  developed  using  the  non-dispersive  asymptote,  with  the  concomitant  assumption 
that  resonant  interactions  are  only  possible  among  coUinear  shallow  water  waves.  This  inherently 
disallows  vector-sum  interactions  among  spectral  components  in  the  wave  field,  negating  a  significant 
portion  of  the  potential  nonlinear  behavior. 

Eldeberky  and  Battjes  (1995)  developed  a  parameterized  triad  energy  exchange  mechanism  which 
depends  on  the  evolution  of  the  biphase,  a  higher-order  statistical  quantity  of  the  wave  train  which 
is  zero  at  low  nonlinearity  (low  Ursell  number)  and  asymptotically  approaches  —n /2  as  the  Ursell 
number  increases.  The  evolution  of  this  quantity  was  determined  from  experimental  measurements 
as  a  function  of  Ursell  number.  Additionally,  nonlinear  interaction  was  further  limited  to  self-self 
interactions  at  the  spectral  peak.  This  allowed  energy  to  move  from  the  peak  to  its  second  harmonic, 
but  did  not  allow  for  any  feedback  transfer.  This  parameterized  model  is  a  selectable  option  in  the 
SWAN  model  (see  Chapter  5). 

Recently,  much  has  been  done  on  developing  stochastic  models  from  the  phase-resolving  dynami¬ 
cal  equations  (for  example,  Boussinesq  equations).  The  evolution  equation  and  its  conjugate  are  each 
multiplied  by  their  corresponding  conjugates,  then  added  together,  resulting  in  an  energy  equation 
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with  triple  products  of  amplitudes  in  the  nonlinear  summations.  These  triple  products  are  related  to 
the  bispectnim  (a  higher-order  spectrum),  for  which  an  additional  evolution  equation  is  needed  for 
system  closure.  The  bispectral  evolution  equation  is  derived  by  applying  the  triple  product  defini¬ 
tion  of  the  bispectrum  to  the  original  evolution  equation,  resulting  in  an  equation  for  bispectra  with 
quadruple  products  of  amplitudes  in  the  nonlinear  terms  (the  trispectrum).  At  this  point  a  Gaussian 
closure  assumption  is  made,  which  has  the  effect  of  reducing  the  trispectrum  to  products  of  the  spec¬ 
trum,  thus  creating  a  finite  system  of  equations.  Herbers  and  Burton  (1997)  developed  a  directional 
stochastic  model  from  the  Boussinesq  equations  of  Peregrine  (1967)  using  the  same  periodicity  as¬ 
sumptions  as  Kirby  (1990)  and  the  procedure  described  above.  Agnon  and  Sheremet  (1997)  worked 
from  the  nonlinear  dispersive  model  of  Agnon  et  al.  (1993)  to  develop  a  stochastic  model  using  the 
same  closure  hypothesis  as  Herbers  and  Burton  (1997),  but  in  the  form  of  a  single  equation  model 
for  the  spectral  energy.  Kofoed-Hansen  and  Rasmussen  (1998)  operated  on  the  extended  Boussinesq 
equations  of  Madsen  and  S0rensen  (1993)  and  developed  a  corresponding  stochastic  model.  They 
showed  that  good  model  comparisons  were  possible  so  long  as  the  Ursell  number  Ur  <  1.5,  was  in 
accord  with  the  observations  of  Agnon  and  Sheremet  (1997).  Eldeberky  and  Madsen  (1999)  revis¬ 
ited  Agnon  and  Sheremet  (1997)  and  rederived  a  stochastic  model,  taking  an  additional  second-order 
effect  (noted  earlier)  into  account. 

One  potential  limitation  has  been  the  specification  of  the  closure  required  to  truncate  the  system 
to  a  finite  number  of  solvable  equations;  the  subset  of  trispectra  used  in  the  bispectral  evolution 
equation  represent  the  lowest  order  contributions.  However,  Holloway  (1980),  in  the  context  of  triad 
interactions  among  internal  waves,  hypothesized  a  different  closure  for  the  system.  Rather  than  dis¬ 
carding  a  significant  portion  of  the  trispectra,  Holloway  (1980)  suggested  that  the  trispectrum  is  also 
proportional  to  the  bispectrum;  this,  in  addition  to  the  products  of  energy  terms,  make  up  the  con¬ 
tributions  from  the  trispectrum.  Becq-Girard  et  al.  (1999),  working  from  the  extended  Boussinesq 
model  of  Madsen  and  S0rensen  (1993),  developed  a  stochastic  model  with  Holloway’s  (1980)  clo¬ 
sure  hypothesis  included.  This  inclusion  essentially  adds  a  linear  term  (multiplied  by  an  empirical 
proportionality  coefficient)  to  the  bispectral  evolution  equation,  with  the  effect  of  broadening  the 
resonance  condition,  and  adding  higher-order  contributions  that  may  improve  overall  performance  at 
moderate  Ursell  numbers. 


8.  CONCLUSIONS 

Nonlinear  frequency  domain  models  have  undergone  rapid  development,  apace  with  correspond¬ 
ing  advances  in  the  time  domain  realm  (particularly  with  respect  to  Boussinesq  models).  They  have 
also  increased  in  utility  with  the  incorporation  of  enhanced  frequency  dispersion  effects,  improved 
shoaling  and  energy  dissipation  from  wave  breaking.  We  investigated  several  formulations  for  nonlin¬ 
ear  frequency  domain  models;  ensuing  data-model  comparisons  demonstrate  that  nonlinear  models 
based  on  the  mild-slope  equations  appear  to  be  more  accurate  than  frequency  domain  transformations 
of  the  extended  Boussinesq  equations  explored  herein. 

Initial  phases  of  the  irregular  wave  train  were  available  to  drive  the  models  for  data-model  com¬ 
parisons  to  the  laboratory  data.  In  most  general  field  situations,  however,  smoothed  spectra  from 
pressure  gauges,  wave  buoys  or  forecast  model  output  would  be  the  only  source  of  data.  Using 
smoothed  spectra  as  an  initial  condition  requires  multiple  runs  of  the  model  with  random  phases,  a 
time  consuming  task.  With  the  advent  of  the  SWAN  model,  the  consideration  of  triad  interactions 
(even  in  the  parameterized  form  used  by  the  model)  has  become  more  widespread,  and  the  need  for 
a  useful  operational  form  of  these  interactions  more  apparent.  This  is  particularly  evident  as  more 
model  systems  linking  wave,  hydrodynamic  and  sediment  modules  are  developed.  The  stochastic 
models  described  previously  exhibit  great  potential  for  operational  use;  more  development  in  this 
area  is  warranted. 
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Z'  —  nonlinear  coefficient  from  Eldeberky  and  Madsen  (1999) 

2"  —  nonlinear  coefficient  from  Kaihatu  (2001) 

z  —  vertical  coordinate 

V  —  horizontal  gradient  operator 
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